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FOREWORD 


This report presents the results of work performed by Computer Sciences 
Corporation while under contract to the Aero-Astrodynamics Laboratory of the 
George C. Marshall Space Flight Center, Contract NAS 8-26113. The program 
described herein evolved as a supplemental application package when developing 
techniques of orbital decay and long-term satellite ephemeris prediction (see 
Ref. 1). 

The author is grateful to Mr. F. C. Boles (formerly with CSC) for his important 
contributions in the initial phase of the program development. 
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SUMMARY 


The symbolic manipulation capabilities of the FORMAC (FORmula MAn ipulation 
Compiler) language are employed to expand and analytically evaluate integrals 
of the form 


I = k 


h 


sm X cos 


Sc 


(1 + e cos x) 


n 


dx 


o < e < 1 


where p, q are non-negative integers and n is an arbitrary integer. Basically, 
the program integration is effected by expanding the integral(s) into a series of 
subintegrals and then substituting a pre-derived and pre-coded solution for that 
particular subintegral. Derivation of the integral solutions necessary for pre- 
coding is included, as is a discussion of the FORMAC system limitations en- 
covintered in the programming effort. 
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SECTION 1 - THEORY 


1. 1 TECHNICAL DESCRIPTION 

Differential equations involving integrands of the type 


k sin^x cos\ 
n 

(1 + e cos x) 


o < e < 1 


(where k^does not depend on x, p and c^are non- negative integers, and n is an 
integer) arise in several areas of general perturbation theory. By properly 
utilizing the algebraic and analytic capabilities of a symbolic manipulation lan- 
guage such as FORMAC, these integrals can be analytically evaluated. The 
program described herein is designed for just that purpose. 

The program consists of a FORMAC driver and a set of subroutines which 
effect the required integrations. The driver performs all required manipula- 
tions of each input integrand, determines the integration parameters £, c[j n 
and the "constant" k, and transmits these quantities to the driver routine of 
the integration package (the set of routines which perform the required integra- 
tions). The integration package driver then identifies the integrand involved, 
makes any necessary variable transformations, and calls upon the proper sub- 
routine to carry out the integration. 


The complete solution of an integrand usually requires solving several sub- 
integrals (special cases), and each integration package subroutine is designed 
to integrate a given type of sub integral. These integrated results are then 
transmitted back to the integration package driver where inverse transforms 
are performed (if necessary), and the results passed on to the FORMAC driver 
for simplification and output. 

A detailed analysis of the required integrations is provided in the appendices, 
and the individual integration subroutines are structured according to this 
development. 
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1.2 EQUATIONS (OMITTING THE INTEGRATION CONSTANT) 
If n < 0 


m . 

I sin^ cos^ (1 + e cos dx = ^ j j sin^c cosSc^ dx 

•' i=0 


( 1 - 1 ) 


where niv= -n and 


f . a b 
sin X cos X dx 


. a+1 b-1 , , 

sin X cos X b-1 

a+b a+b 


a-1 b+1 , 

sin X cos X ^ a-1 

a+b a+b 


/■' 

/ 


. a b-2 , V / „ 
sin X cos X dx b 0 


. a-2 b , . „ ,, , V 

sin X cos X dx a ^ 0 (1-la) 


a + b = 0 


by backward recursion. 

If n > 0 and p is odd, set p = 2K + 1 and x = cos x to obtain 

K 


/ 


. P C 
Sin X cos 


(1 + e cos x) 




K! 


(K - i)! il 


(- 1 ) 


i=0 


d^I 

a (f - 1)! 1£ 

(i + a -1)! ^ a 
de 


( 1 - 2 ) 


where a = 2(K - i) + q < n 
J2 = n - a 


and 


} [l/ e (i - 1)] (1 + ez) ^ ^ ^ £ / 1 

l/e log (1 + ez) j? = 1 


(l-2a) 


(l-2b) 


1-2 



or, 


/ K I b 

sin^x cosSc , V?, i K! )v^ /* b-j 

(1 + e cos x) *' 

' i=0 ( j=l 


dz 


i=0 

n-1 




j=l 

j (£ - 1)! "U 


d^L 


i=o 


D-] ' ' (i + K - 1)! 


where b = a- n+ l a>n 
X = n - j 


1 1 /n + j - 2\ o • u 

ot = - , ot. = - - - . ., — a. ^ 2<]<b 

le ] e\]-l/j-l 


^n-i = -E 


^ n ! e" a 




h-S. 


-j ^ (n - j + J?) ! (K - je) r 0 < j < n - 1 , 0 

je=o 

2 2 

If n > 0 and p is even, set p = 2K and use sin x = 1 - cos x to obtain 


/ sin^ cosSc 
/ (1 + e cos x) 

where a = 2(K - i) + q < n 
J? = n - a 


K 


i = E'>‘ 


K! 


(K - i) ! i ! 


i=0 


(- 1 ) 


.a 

a (jg-- 1)! 2jg 

(Jg + a - 1)! ^ a 
de 


and 




dx 


K(x) + 


(1 + e cos X 


/ 


(1-3) 


(l-3a) 


(1-4) 


(l-4a) 


(l-4b) 


1-3 



with R(x) = 


(1 + e cos X) 


A='{1- 1 ) 


1 - e 


2n - 3 


C = - 


(n-2) 


2 X 

^21= — {(ttt) 

U - e j 


- 1 j e sin X 

(1 - e^) (1 - e^) (1 + e cos x) 


sin^ cos*^ 


(1 + e cos x) 




f H J 
'cos X dx 


Bn-, <-‘> 


i (jg - D! 2jg 
a.K-D! 


where b , J? , a. , oi a . are 

0 1 j n-j 


defined as before. 


1-4 



SECTION 2 - PROGRAMMING 


2.1 GENERAL 

The deck consists of: (1) a main program that determines the integration param- 
eters (p , q , n) for the various integrals appearing in the differential equation(s); 
(2) an integration module driver that executes the proper logic to call the re- 
quired integration routines and accumulates their results; and (3) several 
integration routines that actually perform the required integrations. 

2.2 PROGRAM CHARACTERISTICS 

The program is coded in the FORMAC language and runs on the IBM 7094 under 
IBSYS (Version 13) control. Program input is internal (no card or tape reads) 
and consists of the differential equation(s), along with parameters which indicate 
the initial equation set and how many equations to evaluate per set (provision is 
made for up to six sets of three equations each). Program output (print only) 
consists of the anal 3 rtical integration results and certain parameters which 
provide a trace of the integration flow. 

Program run-time varies from 3 to 20 minutes, depending on the number of 
equations and their individual complexity. Normal output is 0 to 5000 lines and 
is in FORMAC notation (as is the input). One non-system diagnostic appears 
in subroutine MNDVIG. This diagnostic indicates an improperly dimensioned 
array, terminates the current integration, and returns control to the main 
program for continued execution. An A5 work tape is required. 
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2.3 BLOCK DIAGRAM/FLOWCHART 

Block diagrams and flowcharts for the FORMAC main program and its sub- 
routines are presented in Figures 2-1 through 2-6. The blocks depicted in the 
subroutine figures are defined as follows (see Subsection 2.4 for a complete 
description of the subroutines): 


IMl 

IM2 

IMS 

IM4 

T1 

T2 

T3 

R12 

IRl 

IR2 

IRS 


Denotes integration by SPCQ routine. 

Denotes integration by DUIS routine. 

Denotes integration by MNDVIG routine. 

Denotes integration by XN routine. 

Denotes transformation given by 

(sin x)^= (1 - cos^ x)^ ^ ^ ^ 

Denotes transformation given by x = cos x 

Denotes transformation given by x = 1 + e cos x 

Denotes recursion relation 1 or 2 (depending on ITYPE) 
of subroutine DUIS. 

Denotes integration by recursion relation 1 of subroutine 
SPCQ. 

Denotes integration by recursion relation 2 of subroutine 
SPCQ. 

Denotes integration by recursion relation 3 of subroutine 
SPCQ. 
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START 



Figure 2-1, Main Program 















Figure 2-2. Subroutine IDIGTE (Sheet 1 of 2) 
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S2 



S3 



Figure 2-2. Subroutine II 
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(Sheet 2 of 2) 











Figure 2-3. Subroutine MNDVIG 
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SUBROUTINE 

OUlS 

ITYPE,M,N,ANS 



Figure 2-4. Subroutine DUIS 
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Figure 2-6. Subroutine XN 
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2.4 program SUBROUTINES 


All program subroutines are FORMAC routines. 

2.4.1 Subroutine IDIGTE 

This routine uses the integration parameters p , q , n to select the proper 
logic for the type integral being considered and calls the various integration 
routines required to effect the actual integration. 

Argument List Variable Type 

Ml - m (input) FORTRAN 

M2 - n (input) FORTRAN 

ANS - I (output answer) FORMAC 

2.4.2 Subroutine MNDVIG 

This routine performs a multinominal division to decompose the integrand into 
a quotient and remainder and then effects a termwise integration of each. (Two 
basic types of integrals are considered. ) 


Argument List 

Variable Type 

ITYPE 

- selects integral type (input) 

FORTRAN 


values 1 or 2 


M 

- exponent of numerator (input) 

FORTRAN 

N 

- exponent of denominator (input) 

FORTRAN 

ANS 

- integration answer (output) 

FORMAC 


2.4.3 Subroutine PUIS 

This routine integrates by differentiating (parametrically) under the integral 
sign. (Two basic integral types are considered. ) 


Argument List 

Variable Type 

ITYPE 

- selects integral type (input) 

FORTRAN 


values 1 or 2 


M 

- exponent of numerator (input) 

FORTRAN 

N 

- exponent of denominator (input) 

FORTRAN 

ANS 

- integration result (output) 

FORMAC 
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2.4.4 Subroutine SPCQ 


This routine integrates various products of sin x cos x. 


Argument List 

M - exponent of sin x (input) 

N - exponent of cos x (input) 

ANS - integration result (output) 

2.4.5 Subroutine XN 

This routine integrates powers of x . 

Argument List 

N - exponent of x (input) 
ANS - integration result 

2.5 MAIN PROGRAM MNEMONICS 


Symbol 

Mnemonics 

Type 


NTERM 

FORTRAN 


IFLAG 

FORTRAN 

dl/dx 

DEDT 

FORMAC 

(1 + e cos 

X) PBYR 

FORMAC 


SOL 

FORMAC 


SEC 

FORMAC 


NT 

FORTRAN 

X 

X 

FORMAC 


CF 

FORMAC 

P 

P 

FORTRAN 

q 

Q 

FORTRAN 

n 

N 

FORTRAN 


QQ 

FORTRAN 


Variable Type 

FORTRAN 

FORTRAN 

FORMAC 

Variable Type 

FORTRAN 

FORMAC 

Purpose 

Number of differential equations per set 

Specifies initial set of differential 
equations 

Differential equations to be integrated 

Represents expression (1 + e cos x) 

Solution to differential equation 

Secular component of solution 

Number of terms in differential equation 

Integration variable 

Coefficient of integrand 

Exponent of sin x 

Exponent of cos x 

Exponent of (1+e cos x) 

Flag used during expression output to 
signal output termination 
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2.6 PROGRAM LISTING 


The following is a sample of the program listing. 
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SJO^ rll3 MCYFR PTN ^^06»310030.10,12 

SPaUSf RfSfT all rONFOLF KFYS 

5'=:XccUT‘^ IBJOP 

SIRJOP format »F TOFS 

SIRfMT npiVF M94/2 

INTFC-FR output (20) 

integer dummyo 

integer DELDTtE.X , A NS , AN SCK , aS DIJMM yO , DUMMY F 
INTEGfR P ,Q ,DI EXT ,ni iMMY 
intfgfr PDYR 
symaRg 

ATO’^IT M 
ATOMIC PRYR 
ATOMIC F»X»W*r 

LrL Par am (X,M) » (FMCSIN(X ) ,FMcSIN(M)D, (FMccOs(x ) .FMCCCS(M) ) 
LFT ANS=0 
NN = T 

). LfT nUMMYO=F /?» f*fmfS I N ( W ) *FMrrOS ( I ) **? 

L^^T PFLnT =nUMMYO*TM''S ?N(X 1**P*fMcc05(X)*^^^P^YR**? 

LFT DUMMY0=nUMMYO 
GO TO 1001 

lO CALL IDIGTF(2 ,s»?*aNS,LBU 
LFJ ANS=DUMMyO»ANS 
LFT ANS=FXPaND ans 
WRITF(6»'=0) 

WRT TF ( 6 , 1 00 ) 

00 = 0 , 

104 LfT Q0=BCDC0N aNS, output, 19 

WRiTf( 6 » 101 ) (OUTPUT ( KK ) ,KK=? , ion 
IF(OO.NF.0, ) GO TO 1 O 4 

LfT aNS=SUPST aNS , ( fMcS I N ( x ) , FMcS I N ( M ) ) 

LcT aNS=SUPST aNS , ( fMccOS ( X ) , fMccOS (M ) ) 

LfT ANS=SUBST aNS , ( fMtcOS ( M ) »*2 »U 1 -fMCS I N ( M ) ** 2 ) ) 

GO TO 9B 

7 WRTTf( 6 ,fO) 

LfT nUMMYO= ( 1 +f»^<-? ) #* ( 1 / ? ) 

LfT OFLoT=DUMMY 0 *FMrSlN(X)** 2 *FMcc 0 S(X)»* 3 /PRYR **2 
L>^T nUMMYO = DUMMYO 
GO TO ]Oni 

20 CALL IDIGTF(DFLnT,2,3,-2,ANS,DUMMYO) 

GO TO 98 

0 WRTTf(6,fO) 

let nUMMY0=1 

LrT OFLnT = DUMMY0*FMrS IN ( X » »»0 *fMccOS ( X ) **7/PP YR #*0 
LET OUMMYO=DUMMYO 
GO TO 1 001 

30 CALL IDIGTE(DELDT,3,2,-3»ANS,DUMMY0) 

GO TO 98 

4 WRTTf( 6 ,fO) 

LET nUMMYO=l 

|_cT nFLOT = FMCSlN(X)**7/PRYR»»3 
LET nUMMYO=DUMMYO 
. GO TO lOOl. 

4 O CALL IDIGTEIDELDT, 7 , 0,-3, ANS, DUMMYO) 

GO TO 98 
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WRTTcf6,^0) 

LPT nUMMY0=1 

i_ft p.PLnT = PMrsi n ( X ) /ppyr*-^? 

L>^T nUMMY0 = DUMMYO 
GO TO lOm 

'iOl CALL IDIGTE(DELDT,p,0,-2»ANS»DUMMY0) 

GO TO P8 

6 WR T T P ( 6 ♦ A 0 ) 

let nUMMYO=l 

L^T nPLnT=EMrrOS ( X ) **A/PRYR»*A 
LPT nUM'^YO = DUMMYO 
GO TO TOOT 

60 CALL IDIPTE(nPLnT,0,3»-3»AN5»DUMMY0) 

GO TO 98 

98 00=0. 

WRITP(6*80) 

«^0 PORMAT ( 1 HI » / // ) 

'.•'RTTc(6»i 0? ) 

lOp PORMaT ( /// 1 3H COPP’^ICTPNT/ ) 

99 Lt^T 00 = Rr0C0N OUMMYO »OUTPUT * 1 

WRiTf( 6*101) (OUTPUT(II) »II=?»19n 
IP(OO.NE,0.) GO TO 09 

00 = 0 . 

WR T T P ( 6 * 1 0 3 ) 

103 P0RMaT{///3H ANS/) 

100 LFT OQ=RCncON A NS jOUTPUT » 1 9 

WR T Tp ( 6 ♦ I.Ol ) (OUTPUT ( I I) , I I = 2 » 1 on 
IF(OO.NP.O. ) GO TO 1 00 
00 = 0 . 

NN=NN+1 
NN = 7 

GO TO ( 1 * 2 »3 »4 » R *6 »7 ) » NN 
lOOl CONTINUE 

c the pollowing logic is designed to accept a given integrand of THF 

C FORM 

r DTDX = K*STN(X)*»P*cOS(X)**0#(i+pi^cOS(X) ) **N 

C WHERE K IS INDEPENDENT OF X. 

C isolate K.AND DETPRMINF P,Q,N pO- PURPOSES OF INTEGRAL IDPNTIFICAT 
C ION, 

c start logic 

LFT DIDXT=DPLDT 

LFT nuMMY=COFpp D I ny T , PMCS I N ( X ) *M0 , VI » V2 
P=V2 

LfT dUMMY=cOfff D I nxT , fMccOS ( X ) ^*0 , VI » V2 
0 = V2 

LFT DUMMY=cOpfp DIRXT,PrYR»*0»V1,V2 
N = V2 

LFT DUMMY=DIDXT/ ( FMCSI N ( X ) **P*FMCCOS ( X ) **0*PBYR**N ) 

LFT dummy=pxpand dummy 

WRTTF(6*R04^) 

804 FORMAT! 1hO,SHDUMMY,///) 

00 = 0 , 

1026 LPT OQ=RCDCON DUMMY .OUTPUT » 1 9 

WRi T p ( 6»1 01 ) (OUTPUT (KK) .KK = 2 » 190 
1P(OQ,NE.O. ) GO TO l026 ■ 
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WR T TF ( 6 » 1 0«i ) P,0,N 

105 format ( 1h0» IX »2hP = » I 2 » IX ,2hO=» I 2T IX ,2 hN* , f 2/ ) 

C FNH L^GTC 

GO TO ( 10>20, 50, 40,^01,60, 7) ,NN 
lOl format ( IhO, 1946/ ) 

7 STOP 

end 

SI EDIT SYSCKI 

SmFTr riRiVF M94/2 
SIFDIT 

$r«FMC TPMOO MQ4/2 

C SUBROUTINE IDlGTE IS DESIGNED TO INTc^GRA.TE THE FOLLOWING DIFFERENT 
C lAL FORM 

C DI/DX = SlN(X) **P*cOS ( X ) I +F«-COs ( X ) ) 

C P,Q»N integRal,p,o nonnfgativf 

c The solution mfthOd sflectfd depends on p,u,n 

C The INTi^GRaTION is dONf via integration subroutines as R<^QUlREn 

r Thf iNTfgRaTION is dONf via INtfgRaTiON subroutines as ROiilRf-o 

C CALLU'G SEOUEMCF DEFINTTIOA'S 

C INPUT VARIABLE 

c MI ***E-XPONENT OF SIN (X) 

c M?**^^EXP0NFNT OF CQS(X) 

c M5*#*EXP0 NENT of {T+F*C0S(X)) 

C OUTPUT variable 
c ANS***INTEGRAL solution 

c 

SUdROUT INp TDIGTf(M 1 ,M2,M5,aNS,LRL) 

INTEGER A NS, ANSI , ANS I I , ANS I I I , E ,P ,0 ,0 I , X 

SYMARG ANS 

ATOMIC E,X 

AUTSIM QINT 

LET ANS=n 

P=M1 

0 = A ^2 

NsM-^ 

r test n for integrand type 
lE(N.LT.O) go to 2 
c start LOGIC BLOCK FOR N.GF.O 
CALL SPCQ(P,Q,ANS) 

IF(N,EQ.0) GO TO 10 
LET ANSI=0 
DO 1 I =I ,N 
01 = 0+1 
LET ANSII=0 
CALL SPCO(P,OI , ANSI I ) 

LET ANSI=FMCFAG(N)/(FMCEAC( I )*EMcFAC(N-I ) )*F*#I#ANSI I+aNSI 

1 CONTINUE 

let ANS = ANS+ANSI 
tO LET ANS=EXPAN0 ANS 

let ANS=SUBST ANS,LBL 
RETURN 

C END LOGIC BLOCI^ FOR N.GF.O 
C start logic BLOCK FOR N.LT.O 

2 N = -N 

c Particular attention must be paid to all possible cases for this type 
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C INTFGR^^ND IKS THRSF! C.ASFS RFQUIR'^ DIFFFRFNT INTFGRatTON MPTHOHS 
IF( (P.FO.O).ANn. (O.FO.O) ) GO TO F 
IFfP.FO.O) GO TO A 
IF(O.FQ.O) GO TO F 
C DFTFRMINF IF P IS FVFN 

IF( (M0D(P»2) ) .FO.O) GO TO 24 

C SlNCF P IS ODD transform ThF INTFGRAND BY X=rOS(X) TO 
r ni/nX = -M-X** 2 )**K*X**C/Ui+F*X)**N K=(P-n/7 

r each TFRM in the FXPaNsiON OF THF integRano if ORTAINFO, 
c thfn integrated either py difffri fntation under an 
c integral sign or RY multinomial hi VI si on 
K=(P-n /2 

LET aNSI=0 

m6=o 

Ip(MO.GF.N) go to 
CALL OUIS ( 2 »MQ ,N , ANS I ) 

GO TO 22 

71 gall MNdVI G ( 2 *MQ ,N , ANSI) 

22 LET ANS=-ANSI 
IF(K.FQ.O) GO TO 27 2 
LET ANSI=0 

DO 2 7 I = I »IF 
LFT ANSII=0 
MQ=2*I+0 

TP(MO.GF.N) GO TO 770 
call DUIS(2*MQ,N»ANSTI) 

GO TO 771 

230 CALL MNDVIG ( 2 ♦MQ *N» ANSI I ) 

231 LFT ANSI = - ( - 1 ) ** I #FMCF AC ( K ) / ( FMCF AC ( I ) *FMCF AC ( K- 1 ) ) *ANS I I + ANS I 

23 CONTINUE 

LFT ANS=ANS+ANSI 
237 LPT aNS=FXPANd ANS 
LFT aNS=DUMMY»ANS 
LFT aNS=FXPaNd ANS 
return 

24 CONTINUE 

C P IS even, use SlN(X)»*P = ( 1-C0S(X)M#7)**K»K=P/2 
C TO get DI/DX = (1-C0S(X)**7)**K*C0S(X)**0/(1+F*C0S(X) )**N 

K = P/2 
M0 = 0 

LFT ANSI=0 
IFIWO.GF.N) GO TO 7f 
CALL DUIS( 1 »MO»N,ANSI ) 

GO TO 26 

25 CALL MNdVIGI 1 »MQ*N,aNSI ) 

26 LET ANS=ANSI 
IF(K.FQ.n) GO TO 277 
LFT ANSI=0 

DO 7 7 I = 1 , K 

LF.T ansi 1=0 

17=7*1 

MO=l2+0 

IF(MO.GE.N) GO TO 770 
CALL DUI5 ( 1 *MO,N ,ANSI I ) 

GO TO 777 
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P70 tall MNDVIG( ] »MQ*N, ANSI I ) 

271 LET ANS!=(-I ) »# I *FvcFAC ( K ) / ( FMCF AC ( I )#FMCFAC(K-I ) )*ANSri+ANSI 

27 CONTINUE 

LPT aNS=ANS+aN5I 

272 LPT aNS=FXPaND ANS 
LFT ANS=DUMMY*ANS 
LET ANS=FXPAND ANS 
RETURN 

? CONTINUE 

c lNTrr,RANo IS FORM S J N t X ) **P / ( ] +F*rOS ( X ) ) 

c dptfr’^inp if p is pvfn 

IP( IM00(P*2) I .PO.O) GO TO 72 

C SlNfF P IS Odd transform the integrand by X=1+F*c0S(X) 
c TO Dl/nx = (- 1 /P**P)*(P**7-(X-1 )*»2)M^K/x**N 
C WHERE K=(P-1 )/2 

C EACH TFRM in THE EXPANSION OF THE INTEGRAND IS OBTAINEP 

c Them integrated as an imtfgral power of x, 

K= ( P-7 ) /2 
let ansi=o 

CALL XN(-N»ANSI ) 

LFT ANS= (-1/E)*ANSI 

IE(K.EQ.O) go to 710 

LET ANS 1=0 

DO 71 1=1 ,K 

LFT aNSII=0 

17=2*1 

MM=I2-N 

CALL XNIMM.ANSI I ) 

DO 70 J = 1 » I 2 

MMrl 7-J-N 

LET ANSI I 1=0 

CALL XN(MM, ansi I I) ' 

LET ANSI I = ( -1 1 **J*F^<FAC ( I 2 I / ( F^CF AC ( JI^FMCFAC ( I 2~J ) ) * 

*ANSI I I+ANS I I 
70 CONTINUE 

LET aNSI=(- 1 )**( I-l )*F**(-I-I2 )*EMcFAC(<)/( EMCFACI I )*FM rpAC(<-I ' )* 
♦ANSI I+ANSI 
31 CONTINUE 

LET ANS=ANS+ANSI 
710 LET ANS=EXPAND ANS 

LET aNS=DUMMY*ANS 
LFT aNS=EXPaND ANS 

return 

72 CONTINUE 

C P IS PVEN. USE SiN(X)**P=(1-cOS(X)**2)**K*WHFRE 
C K=P/2 to get DI/dX = (1-C0S(X)**2)»*X/( 1+E*C0S(X) )**N 

LET ANSI=0 
K = P/2 

CALL DUIS ( 1 »0»N,ANSI ) 

LFT aNS=ANSI 

lE(K.EQ.O) GO TO 7S2 

LET aNSI=0 

DO 7S 1=7 ».K 

LET ansi 1=0 

17=2*1 
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n ^ 


350 

351 
35 

352 


IP(I2.G'^.N) GO TO ^50 
tall DUIS(1»I2»N»ANSII) 

GO TO 351 

CALL MNDVIG ( 1 ♦ I 2 »N » ANSI I ) 

LFT ANSI = (-1 ) #*I»FMCFAC(K ) /( FMCFAC( I ) *FMCFAC( K-I ) ) »ANSt I + AN SI 

CONTINUE 

let ANS=ANS+ANSI 

lft ans=fxpand ANS 

LPT ANS = nUMMY<fANS 
LPT aNS=FXPaND ANS 


return 

CONTINUE 

INTfGRaNO is FORM ni/DX = C I S < X ) **Q/U 1 +F»cOS ( X ) ) *N 
LET ANSI=0 
IF(Q.GP.N) go to AO 
CALL nuis ( 1 ♦0*N,ANSI ) 


GO TO 41 


40 gall mnovig( 1 »o»n» ansi ) 

41 LET ANS=ANSI 

let aNS=EXPaND ANS 
let ans=dummy»ans 

LET ANS=EXPaND ANS 
RPTURN 

5 CONTINUE 

C integrand is form DT/DX = 1 /(1+F»C0S(X) )**N 

let ansi=o 

CALL nuis ( 1 .0,N.aN5T ) 

LFT ANS=EXPaNO ANSI 

let ans=dummy*ans 
Let ans=expand ans 

RETURN 

C END LOGIC FOR N.LT.O 


END 

SIEDIT SYSCKI 

$I«3FTC IPM06 M94/2 
SIEDIT 

SIBFMC IPMOl M94/2 

C SUBROUTINE MNQVIG INTEGRATES TWO DIFFERENT INTEGRaND 
C FORMS BY multinomial DIVISION* FOLLOWED BY TERM 
C WISE INTEGRATION OF A FINITE SUM OF BASIC INTEGRAND 

c the integrand forms ARF 

c DI/DX=-rOS(X)**M/( l,+P*cOS(Xn*MN M.GF.N 

c DI/DX = X**M/ { 1 ,+F*X )**N M.GFCN 

c the division yields a quotient wh+ch is a 
C SUM OF (M-N+1) basic INTEGRANDS OF T+E FORM 

C A( J)*C0S(X)»*(M-N+1-J) 

C OR A( J)*XiHf (M-N + 1-J) 

c with The a(Ji constant aNd 1 , Lf, J . LFC < M-N+1 ) * 

C AND A REMAINDER WHICH IS A SUM OF N BASIC 
c integrands of The form 

c B(N-K)*C0S(X)**K/( 1 ,+P*COS(X) )*MN 

C OR B(N-IC)*X**K/(1,+E*X)**N 

c with The b(n-ki constant and o.le.k.le. (N-i ) 

C CALLING sequence DEFINITIONS 
C INPUT VARIABLES 
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c 

c 

ITYPf 

*** 

INDICATES integrand TYPf. 
ITYPF.FQ.1 for type 1. 2 FOR TYPF2 

c 

M 


FXPONENT OF NUMERATOR 

c 

N 


EXPONENT OF DENOMINATOR 

c 

ANS 

** 

integral SOLUTION 


c 

SURROUTiNF MNDVIG( ITYPE*M»N,aNS) 
integer A» ANS» ANSI , ANSI I »R .RSUM ,f»X 

symarg ans 

ATOMIC R»X 

dimension A( ?0) *R(?0) 

autsim qint 

LET ANS=0 
MMrM 
NN = N 

MAXA=MM-NN+1 

ITYP=ITYPE 

C ARE A, 8 COEFFICIENT ARRAYS LARGE ENOUGH 
IF(MaXA.GT,;>0) go TO 6 
IF(NN.GT,20) go TO 7 
C SET UP A COEFEICIENT ARRAY 
LET A(1 )=1/e**NN 
IF(MaXA,EQ.1 ) GO TO 10 
00 1 I=2,MAXA 

LPT A(I) = (-l/P)*(NN+I-2)/(I-n*AUI-l ) 

1 CONTINUE 
10 CONTINUE 

C set' UP B COEFFICIENT ARRAY 
let P(NN)=-A(MAXA) 

Nl=NN-l 
DO 2 K=T»NI 

LFT RSUM=FMcFAr ( NN ) / 1 PMCF A r ( NN-Ko*FMCFAC ( < ) ) *F»*K*A( M axA ) 

DO 7 L=T»K 

IFIMaXA.LE.LI GO TO 3 

let R5UM = RSUM + FMCEA.r (NN ) / ( FMCP ACUNN-K + L ) *FMCFAC (<~L ) ) K-L ) *A f M 

*AXA-L) 

3 CONTINUE 

NNK=NN-K 

let B(NNK )=-8SUM 
LFT BSUMrO 

2 CONTINUE 

c integraTf quotient termwise 

let aNSIsO 
DO A 1 = 1. MAX A 

maxai=maxa-i 

IF( ITYP.EQ.2) GO TO 40 
call spcQio.maxai .ansi ) 

GO TO 41 

40 CALL XNIMAXAI .ANSI ) 

41 LFT aNS=aNS+A(I )*ANSI 
LET A.NSlrO 

4. CONTINUE- 

c INTEGRATE Remainder termwise 
let ansi=o 

CALL DUISUTYP.O.NN.ANSII 
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LFT ANS=ANS+B (NN ) *ANST 
IF(NI.FO.O) return 
DO 5 K=l»Nl 
let ANSI=0 
KI =K 

CALL DUIS( iTYPjKl .NN.ANSI ) 

NNK=NN-K 

LFT aNS=ANS+R ( NNK ) #aNSI 

5 CONTINUE 
RETURN 

6 WRtTf(6*60) 

60 format (/// 17H A aRRaY TO SMALL) 

RETURN 

7 WR I T E ( 6 » 7 0 ) 

70 format (///17H B array TO SMALL) 

RETURN 

END 

SirniT SYSCK1 

SIRFTC IPMOl M94/7 
siedit 

SI3FMC IPM02 M94/2 

c subroutine duis integrates The following forms 

c DI/DX= cOS(X)*#M/(t+E*COS(X) )**N M.LT.N 

c nT/OX= X<^*M/ ( 1 +E^^x ) **N M.LT.N 

C The INTEGRAL EVALUATIONS REQUIRf ThE USE Op RECURSIVE 
C RELATIONS ANd/OR D T FFFRENT I AT+ON UNofR AN INTEGRAL 

C SIGN. 

C CALLING SEQUENCF DEFINITIONS 

C INPUT VARIABLES 

C ITYPE *** selects INTFGRaND TYPE 
C M *** EXPONENT Op NUMERATOR 

C N *** EXPONENT OF DENOMINATOR 

C OUTPUT variable 
C ANS *** -INTEGRAL SOLUTION 

C 

SUBROUTINE DUIS( TTyPE»M,N»aNS) 

integer a ♦ ARG *ANS11 * aNS12 »ANS11 I TANS12I * aNS 1 3 . aNS2 . R , C » F , X 

integer ANS 

SYMARG ANS 

ATOMIC E»X 

AUTSIM QINT 

LPT ANS=0 

LFT ANSI 11=0 

let ANS12I=0 

let ANS2=0 

MM = M 

NN = N 

ityp=itype 

K=NN-MM 

c dftfrw.tnf integrand type 
IF( ITYP.EO.2 ) GO TO 2 
C Set up basic SOLUTIONS pOR TYPp 1 

let ARG=( ( 1-F)/(1+F) )*«■( 1/2)»FMcSIN(X/2)/FMccOS(X/2) 

LFT aNS 11 = 2/ (1-F**2 )**n/2 )»FMCATN( aRG) 

LfT ANSl2=l/t 1 -F**P ) *aNSI l-P*FMrSlN ( X ) / ( ( 1 ~F**2 ) * ( 1 +F*FMrcOS ( X ) ) ) 
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IF(K.GT,?) go to 1? 

IFIK.FO,?) GO TO TO 
1F(MM.6T,0) GO TO 1 
LFT aNS=ANS11 
RFTURN 

1 LFT aNS=(- 1 )*»MM*FMrFAC(K-l ) / fMCFAC ( K+MM-1 ) *FMcd I F ( ANS 1 1 ,F*MM) 
RFTURN 

10 IF(MM,GT,0) GO TO 11 
LFT aNS=ANS1? 

return 

11 let AN5= (-1 )**MM»FM.CFAC( K-1 ) /FMCFAG ( K + MM-1 ) *FMcn j F ( ANS 1 2 » F ,MM ) 
RFTURN 

12 LET A = -(NN-1)*(],-E**2)/F 
let b=(2*NN-A)/E 

LET C=-(NN-7)/E 
LET ANSUI=ANS11 
LFT ANST2 I=ANS1 2 
LFT ANSTA=0 
DO 12 1=2. K 

LFT ANSl2=ANS12+FMrSlN(X ) / ( 1 +F*FMCCOS ( X ) )**(NN-1 ) 
»-r*ANS12I-C*ANS11I 
LET ansi, 2 = ANSI 3 /A 
LFT ANS11I=ANS12I 
LET aNS12I=ANS12 
12 . GONTTNUF 

IF{MM,EO.O) GO TO 1A 

let aNS=(- 1) **MM*FMrFAG f X-1 i /r Mf FAG < K+MM-1 ) #pMCniF ( ANS13 >MM ) 
RETURN 

14 LET aNS=ANS12 
RETURN 

2 CONTINUE . 

C CHECK FOR K.FQ.l 

IFIK.EQ.l ) GO TO 21 

LcT aNS2=1/( (T-K)*c)#(i+f*X)**( 1 -K) 

IFfMM.EO.O) GO TO 20 

Let aNS»(-1) **MM*FMCFAC ( K-1 ) /FMCFAC < K+MM-1 ) »FrMcniF ( ANS2 »F ,MM 1 
LFT aNS=SUBST aNS. ( X »FMCC0S( X ) ) 

RETURN 

20 let ANS=ANS2 

LFT ANSrSUPST ANS. (X.FMCCOS(X) ) 

RETURN 

21 LFT ANS2=1 /f*FMCLOG( 1+F*X ) 

' IFfMM.EO.O) GO TO- 22 

let ,aNS=(- 1 )**MM + FMrFAC(K-l)/FMCFAC(K+MM-l )*FMrDIF( ANS2 .F.MM) 
LFT aNS=SUBST ANS. (X.FMCCOSfX) ) 

RETURN 

22 LFT ANS=ANS? 

LFT aNS=SUBST ANS. f X .FMCCOSf X ) ) 

return 

END 

SIEDIT SYSCKI 

SIPFTC IPM02 M94/2 
SIEDIT 

SIRFMC 1PM02 M94/2 

C subroutine SPcO integrates THE FOLLOWING FORM 
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r HT/nx = s I N ( X ) **M*rOs ( X ) M»N,r:F,o (iNtfgfR) 

C STANDARD RECURSIVE RELATIONS ARF USED IN THE EVALUATION 
C CALLING SEQUENCE DEFINITIONS 
C INPUT VARIABLE 

C M *** EXPONENT Op 5IN(X) 

C N *** exponent OF COS(X) 

C OUTPUT VARIABLE 

C ANS *** integral solution 

C 

SlJqROUTlNc SPcQ(M.N,aNS) 

integer a * ANS » ansi tF » RSAVE »X 

SYMARG ANS 

ATOMIC X 

MMsM 

NN = N 

LET ANSsO 

TF{ (MM.fO.O) .and. (NN.FQ.O) ) GO TO A 
IF(MM.FO.O) GO TO 7 
IF(NN.EO.O) GO TO 7 
LET R=1 

I Let ans=ans+b/(Mm+nn)*fmcsin(X)*m(mm+i )*fmccos(X)**(Nn 

LET A= (NN-1. ) / CMM+NN ) 

LPT R=B*A 
NN = NN-;? 

IFfNN.GT.O) GO TO 1 
IF(NN.LT.O) RETURN 
LET PSAVE=B 
LET R=T 
LET ANSIsO 

lO LET aNSI=ANSI-R/MM*eMcSIN(X)**(MM-1 )*FMCC0S(X) 

LET A= (MM-1 ) / (MM ) 

LET B=B*A 

IF(MM.GT.O) go to 10 
IF(MM.LT.O) GO TO 1 1 
LET ANSI=ANSI+B*X 

II LET ANS=ANS+RSAVE*ANSI 

RETURN 

? let P=1 

?0 LET ANS = ANS+R/NN*FMrSTN(X)*FMrcOS{X)**(NN-l ) 

LET A=(NN-1 )/(NN) 

LET P=R*A 
MN=NN-2 

lE(NN.GT.O) GO TO ?0 
lE(NN.LT.O) RETURN 
LET ANS=ANS+B*X 
RETURN 

7 LET R=1 

70 LET ANS=ANS-R/MM*EMrSTN(X)**(MM-l )»FMcrOS(X) 

LET A=(MM-1)/MM 
LET P=B*A 
M.M=MM-2 

IF(MM.GT.O) GO TO 70 
IF(MM.LT.O) RETURN 
LET ANS=ANS+B*X 
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RFTURN 

4 L^T aNS=ANS+X 

RFTURN 
FND 

SlFniT SYSCKi 

SIRFTC IPMOR M94/2 
SIEOTT 

SiRFMr TPM04 MQ4/? 

c SUrROUTINf XN iNTfgRMcs all intfgral POWFRS Of X 
C CALI.IMO sfqufncf DFFINTTIONS 
C INPUT VARIABLE 
C N*** EXPONENT OF X 

C OUTPUT variable 

C ANS*** INTEGRAL SOLUTION 

C 

subroutine XN(N,ANS) 
integer ANS,X»F 
symarg ans 
atomic x»e 

LET aNS=0 
NN = N 

IF ( (NN + l ) ,EQ,0) GO TO 1 
NPI=NN+1 

LFT ANS=X**NPi /NP i 

LfT aNS=SUrST aNS. ( X » n+F*rMrcOSUX ) ) ) 

RFTURN 

I LFT ANS=FMCLOG (X ) 

let ans=subst ans, ( X . ( 1 +E*FMCC0S( X ) ) ) 

RFTURN 

END 

SI ED IT SYSCKi 

SIBFTC IPMOA M9ii/P 
$ IFDIT 
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2.7 DECK SETUP 


Figure 2-7 represents the deck setup of the program and Figure 2-8 is a 
sample instruction card. 


J 
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Figure 2-7. Deck (Sheet 1 of 2) 
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f 7/8 
I'^SIEDIT 
f SIBFTC IPM04 


SIEDIT 


/ 

f SUBROUTINE XN 

[ SIBFMC IPIVI04 

^ SIEDIT 

SIBFTC IPM03 



Figure 2-7. Deck (Sheet 2 of 2) 
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IF EXCEEDS MAX: 


DSTR aSTi QOMPQRETSYI 


FAST TAPES ! A B C 0 


INPUT TAPES 


LOGIC REEL NO. | OEN 



LINES OF OUTPUT (lOOO'S) time. 

-5 a 5- 15 aiS-30 QOVER HOURS 


>ROGRAMMER COMMENTS: InuMBER OF CASES 




CrypicnLV 


OPERATOR COMMENTS: 


a SEE ON-LINE 
a SEE TECHNIQUE! 
□ MAX EXCEEDED 
I D RETURN TO SYS 
D LINE MAX 


OPEN INIT;. 

over;- 


OUTPUT TAPES ONLY 



UNIT NO OF CPYS SAVE 



NO FILES I NO FRAMES COPIES DENSITY COPY-FLO p^^YAR 


MSFC - Form 533 (Rev February 1966^ 


Figure 2-8. Sample Instruction Card 
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2. 8 LIMITATIONS ENCOUNTERED 


The application of FORMAC, as described herein, revealed several system 
limitations. The first, and most severe, of these limitations was the lack of 
sufficient free list storage. On a number of occasions, program execution was 
terminated while expression manipulation was being carried out. This required 
breaking down sequences of analytical computations to obtain partial solutions 
and, on some occasions, nmning the program "piecemeal." 

Problems were also encoimtered in communicating between the various program 
elements, i. e. , between the main program and a subroutine, or between two 
subroutines. For example, a result obtained in a subroutine could be communi- 
cated to the main program for I/O purposes, but could not be differentiated. 

The isolation of the integration parameters p , q^ , _n was unduly complicated 
since there was no convenient way to separate the numerator of a fraction from 
its denominator. In conjunction with this, it was discovered that the automatic 
simplification transformations did not work as expected. For example, one 
such transformation is stated as 


log (X . y) = log X + log y 


Since expressions are stored in internal Polish delimeter form (Ref. 2), it was 
expected that application of the transformation to a product expression such as 

r 3 2I 

log Isin x/(l + e cos x) J 

would yield 


3 log sin X - 2 log (1 + e cos x) 


However, it did not; in fact, no transformation was made at all. The only 
explanation for this is that must be separate expressions before the 

transformation can be applied. 
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These limitations, while not entirely prohibitive, were an inconvenience and 
were circumvented only by special programming efforts. Such efforts required 
additional programming time and detracted slightly from the overall program 
efficiency. 
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APPENDIX A - INTEGRAL SOLUTION FORMULATION 


This appendix contains the formulation of the solution for integrals of the form 


/: 


sin^ cosSc dx 


(A-1) 


(1 + e cos x) 

where £ , q , £ are integers (p , q > 0) and 0 < e < 1 . Two basic cases 
arise, depending on whether n_ is negative or non-negative. 


Case 1 (n negative). If n is negative the integral becomes 

QinBc cos Sc (1 + e cos x) ^ dx 
and application of the binomial expansion transforms this integral to 




(A-2) 


m 

E m! e^ / • P q+i j 

= — : — r I sm X cos X dx 

(m - i)! i!y 


m = -n 


( 1 - 1 ) 


i=0 


Each integral appearing in this summation can be evaluated by repeated applica- 
tion (backward recursion) of 


/ 


.a b , 
sin X cos X dx = 


sinSc ^ cosHc ^ b - (2K - 1) / . a b-2K ^ 

^ ; I Sin X cos X dx 

a + b a + b 


fs 


(l-la) 


K = 1, 2, ... , until b - 2K = 1 or 0 . 

Case 2 (n non- negative). If £ is non-negative, two subcases arise as £ is 
even or odd. If £ is odd, transform the integrand by z = cos x and employ 
the binomial expansion to yield 

-E->‘ 

i=0 


^^2(K-i)+q 
7 (1 + ez)' 


■dz 


(p = 2K + 1) 


(A-3) 


A-1 



If p is even, the binomial expansion (using cos x = 1 - sin x ) will yield 


K 



i=0 


(K 


K ! / cos 

^ 7(77 


2^K-i)+q 


(1 + e cos x) 


dx 


(P = 2K) 


(A-4) 


In either case, the required integrals can be obtained (see Appendix B) by 
parametrically differentiating (Ref. 3) one or the other of two fundamental 
integrals, or by a multinomial division followed by a termwise integration of 
the quotient (integrating powers of z or cos x ) and a termwise integration of 
the remainder by parametric differentiation of these same fundamental integrals. 
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APPENDIX B - DERIVATION OF THE FUNDAMENTAL INTEGRALS 


This appendix contains the derivations of the solutions for the two integrals 

(z = cos X) 


and 






a j 
z dz 


(1 + ez) 


n 


'2 = 


a j 

cos X dx 


(1 + e cos x) 


n 


required in Case 2 of Appendix A. The solution procedures are essentially the 
same in both instances. 

Solution for I ^ . If a < n , define the fundamental integral I^^ by 


■u = 


h 


dz 


Ji = n - a 


(1 + ezr 

and differentiate 1^^^ ^ times with respect to the parameter to obtain 


= (-1) 


g - 1)1 


L(Ji + a - 1)1 


d^I 


11 


de 


Ij^jj is easily obtained as 

1 


e(i - 1) 


(1 + ez) t^l 


^Ul = 


f— log (1 + ez) 


£ = 1 


If a > n , expand the denominator and divide to obtain 


a 


a-n+1 n-1 

a-n+l-j 


(1 + ez) 


a-n+l-j 


K 


j=l 


K=0 


(1 + ez) 


(B-1) 


(B-2) 


(l-2b) 


(B-3) 


B-1 



where 


„ = 2<j<a-n+l 

1 e ’ j e \ j - 1 / j-1 


^n-K=- Z 


K , K-je 
n ! e at 


a-n+l-J? 


(K - £)! (n - K + Jf)!] 


0<K<n-l,« =0 
o 




Then 


a-n+1 n-1 

“ a-n+1- j 


‘i = Z“j 


I, = > a. I Z ^ dz + 

3=1 

a-n+1 


J. jm 

2>a-K 

K=0 J 


= V “i 

a - n + 2 - j 


a-n+2-j 


j=l 


Ba-K 

K=0 

(J? - 1)! ‘^^l 

(I.K-l). 


(l-3a) 


(B-4) 


where j? = n - K and is given above. 

Solution for I - . The solution procedure here is entirely analagous to that 
employed for 1^^ , but slightly more complicated. If a < n define the funda- 
mental integral 


^2je. 




dx 


j 2 = n - a 


(B-5) 


(1 + e cos x) 


and differentiate ^ times with respect to the parameter ^ to obtain 


I2 = (-1) 


a - 1 )! 


(je + a - 1)! 


d^I 


2j? 


(B-6) 


de 


B-2 



B-3 


lo lx 



In de- 



(l-4c) 


(B-7) 


(l-4b) 


(l-4c) 



a-"”--) 


SO that I is known for all 


If a > n 


a-n+1 


E",/ 


a-n+1- j , 
cos X dx + 


n-1 r 

S-k/ 

K=0 J ' 


cos X dx 


(1 + e cos x) 


a-n+1 


s,/ 


a-n+l-j , 
cos X dx + 


I>.-K 


K (Jg - 1)! 2 jg 

(jg + K - 1)! , K 
L de 


where a. , B ,, are defined as before and 
] n-K 

/* m , 1 m-K m-K f m-2K , 

/cos X dx = — cos X Sin X + /cos x dx 

J m m J 


K = 1, 2, • • • , by backward recursion until m - 2K =1 or 0 . 
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